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ABSTRACT 

We outline the steps needed in order to incorporate the evolution of single and binary 
stars into a particular Monte Carlo code for the dynamical evolution of a star clus- 
ter. We calibrate the results against TV-body simulations, and present models for the 
evolution of the old open cluster M67 (which has been studied thoroughly in the liter- 
ature with iV-body techniques). The calibration is done by choosing appropriate free 
code parameters. We describe in particular the evolution of the binary, white dwarf 
and blue straggler populations, though not all channels for blue straggler formation 
are represented yet in our simulations. Calibrated Monte Carlo runs show good agree- 
ment with results of iV-body simulations not only for global cluster parameters, but 
also for e.g. binary fraction, luminosity function and surface brightness. Comparison 
of Monte Carlo simulations with observational data for M67 shows that is possible 
to get reasonably good agreement between them. Unfortunately, because of the large 
statistical fluctuations of the numerical data and uncertainties in the observational 
data the inferred conclusions about the cluster initial conditions are not firm. 
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1 INTRODUCTION 

The mod elling of individual globu lar clusters has a long his- 
tory (see iMevlan fc Heggie! (|l997l ). especially §§3 and 7.7). 
Much of the focus of this work is on static models such as 
the King model and its variants. In this kind of modelling 
the dynamical history of the cluster is almost irrelevant, ex- 
cept for the general assumption that the cluster is almost 
relaxed. By contrast, there have been a small number of 
studies based on techniques which can follow the dynamical 
evolution of a cluster. Most of this work has been performed 
with a Fokker-Planck scheme using finite differences, but 
also there are examples of the use of fluid and Monte Carlo 
methods (Tab.[TJ). 

In the present paper we develop the Monte Carlo tech- 
nique further, and apply it to a new object. The dynami- 
cal ingredients of the Mont e Carlo code are essentially the 
same as those described in iGiersd (120061), whose code em- 
bodies several features introduced by Stodolkiewicd {l986), 



wh ose code was in turn based on that originally devised 
bv lHenonl i|l97llh Three main features distinguish the code 
which is described in t he present paper from that used by 
iGiersz fc Heggie] (I2003T ) in their work on to Cen: (i) it now 
incorporates dynamical interactions between binary and sin- 
gle stars, between pairs of binaries, and interactions of three 
single stars resulting in the creation of new binaries, all 
using cross sections; (ii) it re places the skeletal approach 
to stellar evolution tak en fromlChernoff fc Weinberg (Il990l ) 
by the algorithms of IHurlev. Pols fc Tout! (|2000l ) for the 
evolution of single s t ars, s upplemented by the methods of 
IHurlev. Tout fc Polj <|2002f ) for the internal evolution of bi- 
nary stars; (iii) a better treatment of the escape process in 
the presence of a static t idal f ield according to the theory 
proposed by iBaumgardtj (|200lT ) . 



There are several factors which motivate this work. Star 
clusters are the focus of severa l intensive observational cam- 
paigns (e.g. iBedin et al.ll200ll: iBedin. Piotto fc Kingl | 2003 |; 



Grindlav et al.ll200ll; iPiotto et al.ll2002l; iKaliari et al 
Kafka ct al. 2004; Richer ct al. 2004: Anderson ct al. 



2003 



20061 ). 
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which are now turning to an examination of the parameters 
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Table 1. Dynamical evolutionary models of individual globular clusters 



Cluster 


Method 


Reference 


47 Tuc 


Fokker-Planck 


Behler et al. (2003). Murohv et al. (1998) 


M15 


Fokkcr-Planck 


MurDhv et al. (2003. 1997. 19941. Dull et al. (19971. Grabhorn et al. (19921 


M30 


Fokker-Planck 


Howell, Guhathakurta & Tan (20001 


N6397 


Fokker-Planck 


Dull et al. (19941. Drukier (1993. 19921 


N6624 


Fokkcr-Planck 


Grabhorn et al. (19921 


ui Cen 


Monte Carlo 


Giersz & Heesrie (20031 


M3 


Moment equations 


Angelctti, Dolcctta & Giannone (1980) 



of their populations of binaries and blue stragglers (BS). 
Dynamical models are needed for the design and interpre- 
tation of observational programmes: how is the period dis- 
tribution and the spatial distribution of binaries affected 
by dynamical evolution? Another problem is the abundance 
and spatial distribution of blue stragglers, which can only 
be answered by a technique which follows simultaneously 
both their dynamics and internal evolution. While A-body 
techniques may ultimately be the method of choice for such 
studies, systems with the size of a globular cluster are likely 
to remain beyond reach for some years, simply because of 
the number of stars and the population of binaries. After all, 
it is only recently that the "hardest" open clusters have been 
modelled at the necessary level of sophisticati on, and even 
then the typical simulation takes one month l|Hurlev et al.l 
l2005h . These authors focused on the old open cluster M67, 
which has been chosen by the MODEST intern ational col- 
labor ation (Modelling DEnse STellar systems) (Sill s et all 
120031 ) as a target cluster for comparison between observa- 
tions and various techniques of numerical simulation. We 
also focus on this cluster, partly for the purpose of refining 
our calibration of the Monte Carlo method. 

This paper begins in Sec. 2 with a summary of the fea- 
tures which have been added to the Monte Carlo scheme. 
We also show there how we calibrate the Monte Carlo tech- 
nique with A^-body simulations. Next (Sec. 3) we apply the 
technique to construct a dynamical evolutionary model of 
the old open cluster M67, and compare our results with ob- 
servations. We give predictions for the initial parameters of 
the old open cluster M67. The final section summarises our 
conclusions, and discusses some of the main limitations of 
our models. 



2 TECHNIQUE 

2.1 Coding of binary- and single-star evolution 

From the dynamical point of v iew our Mont e Carlo code is 
almost exactly as described in iGiersj (|2006l ) . In this tech- 
nique a star cluster is treated as a collection of spherical 
shells, each one representing a single star with a certain en- 
ergy and angular momentum. Neighbouring shells are al- 
lowed to interact and exchange energy and angular momen- 
tum at a rate determined by the theory of relaxation. Es- 
capers are removed according to a prescription which mim- 
ics the effect of a tide. Shells corresponding to binary stars 
also interact with single stars, and other binary stars, at 
rates determined by cross sections drawn from the litera- 
ture |Giersj|200ll ). The only dynamical alterations deal with 
tightly bound subsystems, which often arise in systems with 



a large mass range (e.g. those including both stellar-mass 
black holes and stars at the hydrogen-burning limit of the 
main sequence). 

The introduction of stellar and binary evolution has 
been greatly facil itated by the development of the "Mc - 
Scatter" interface jHeggie. Portegies Zwart fc Hurlevl2006h . 
which provides subroutines for initialising the stellar evo- 
lution of single and binary stars, and for retrieving the 
results of subsequent evolution, mass loss, merging of bi- 
nary components, etc. At present two such packages for 
stellar evolution can be employed . One of these is SeBa 
|Portegies Zwart fc Verbunj 1 19961 ) , which is i ncorporated 
within the STARLAB environment l|Hutll2003l ). The other 
is referred to as "BSE" (binary star evolution), and is based 
on the extensive formulae for the evo lution of single stars 
of a range of metallicities given by iHurlev. Pols fc Tout! 
(2000 ), along with the tre atment of binaries presented by 
IHurlev. Tout fc Pols! (|2002l ). Most of our effort has been con- 
ducted with BSE, partly to minimise any development prob- 
lems with mixed-language programming, and partly because 
SeBa is at present restricted to solar metallicity, whereas 
our interest is mainly directed to globular clusters. Gener- 
ally speaking, the use of the McScatter interface poses few 
problems: 

(i) There was one instance of a named common block in 
BSE which by coincidence was the same as the name of one 
common block in the Monte Carlo code; a change of name 
in the Monte Carlo code was sufficient cure. 

(ii) The enumeration of stars requires care. Although it is 
not clear from the interface, the numerical identity of each 
binary determines the numerical identity of the two single 
stars from which it is composed, and it is important that 
the numerical identities of all single stars (both binary com- 
ponents and those which are genuinely single) are different. 

During a time step of the Monte Carlo code, the changes 
caused by relaxation and dynamical interactions between 
binaries and single stars are performed, and then the stellar 
evolution of all stars and binaries is updated. The associated 
loss of mass (if any) is incorporated into the data for each 
star and binary in the Monte Carlo code, and any mergers 
are dealt with by altering the numbers of single and binary 
stars and adjusting the parameters of the bodies affected. 

2.2 Calibration 

In a Monte Carlo simulation it is usual to adopt units 
such that the constant of gravitation, the initial total mass 
and the initial virial radius are 1. In order to incorpo- 
rate stellar evolution into a Monte Carlo simulation, di- 
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Figure 1. Comparison of TV-body and Monte Carlo simulations 
for models with tidal cutoff. The mass (normalised by its initial 
value) is given as a function of time. The initial conditions are 
given in Tab. [2] and described in the figures. The value of TV 
referred to is N s + TV;,, but this differs from the total number of 
particles ( N s + 27V;, ) . The TV-body model is the heavy continuous 
line, and the others are Monte-Carlo simulations with (from the 
left) 7 = 0.04, 0.03, 0.02, 0.01. 



a 1 




Figure 2. Comparison of the evolution of the half- mass radius in 
TV-body and Monte Carlo simulations for models with tidal cutoff. 
The initial conditions are given in Tab.|2]and described in the fig- 
ures. The value of TV referred to is N s + Nt,. The TV-body model is 
the heavy continuous line, and the others are Monte-Carlo simu- 
lations with (from the left, at late times) 7 = 0.04, 0.03, 0.02, 0.01. 
The mismatch at early times in the upper panel is discussed in 
[2~2~2l 



mensional values for the initial total mass and virial ra- 
dius (or equivalent) must be specified, and then the unit of 
time in the Monte Carlo code (which is essentially a relax- 
ation time) is expressible dimensionally with a factor pro- 
portional to 7V/log(7TV), where 7 is a constant and TV is 
the number of stars in the system. While the value of 7 
is rather well known for the case of single stars of equal 
mass, i.e. 7 = 0.11 approximately (iGiersz fc Heggid 1 19941 : 
Ijoshi. Rasio fc Portegies Zwartl |2000| ). the case of unequal 
masses with a population of primordial binaries has been 
studied much less. On analytical grounds Irlenonl l)l97ot ) gave 
a formula which, by way of example, yields a value 7 = 0.007 
approximately for a mass function for single stars of the form 
dN oc m~ 2 dm, over a range in which the mass ratio between 
the maximum and minimum mass is 500 : 1. This theory also 
implies that the value depends o n the mass ratio, which 
changes through stellar evolution. iGiersz fc Heggid (|l996l ) 
found a value 7 ~ 0.015 for a power-law mass function of 
index —2.5 and a smaller mass range of 37.5 : 1, by means 
of intercomparison of TV-body simulations, and somewhat 
larger values from Henon's formula or from comparison with 
isotropic Fokker-Planck models. 

Here we adopt a pragmatic approach, comparing Monte 



Carlo and TV-body models with identical initial conditions. 
These are summarised in Tab.0 where the tidal radius refers 
to a tidal cutoff (or, for the TV-body models discussed from 
S ec !2. 2. 21 onwards, a tidal field). It is necessary to carry out 
this comparison for at least two values of TV. If we were 
to determine 7 from a single value of TV, it might be that 
this value simply obscures some systematic problem with 
the Monte Carlo code, and would fail for a different value 
of TV. Modest values of TV are better for this purpose, as the 
TV-dependence of the Coulomb logarithm becomes weaker 
as TV increases. We study cases with TV = TV, + 27V(, = 3750 
and 15000, where TV s ,TV b are the initial numbers of single 
and binary stars, respectively. 

The Monte Carlo code free parameters are as foll ows: (i) 
7, (ii ) fimin , minimum value of the deflection angle ijCiersd 
1998), (iii) r, the overall time step and (iv) a, see for the 
definition Sec. f/X2~2l 

2.2.1 Models with tidal cutoff 

First, we concentrated on calibration of Monte Carlo models 
for which the influence of the tidal field of a parent galaxy is 
characterised by the tidal energy cutoff - all stars which have 
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Table 2. Initial conditions of calibration runs 



N 3 + 2N b 
Initial model 
Initial tidal radius 
Initial half-mass radius 
Initial mass function 



Binary fraction 
Binary eccentricities 
Binary semi-major axes 



3750 (15000) 
Plummer 
30pc (30pc) 
3pc (3pc) 
Kroupa, Tout . 



Gilmore 



with ai = 1.3, mass range 
between 0.1M Q and 50Af Q 
N S /(N S +N b ) = 0.5 
/(e) = 2e 

Uniformly distributed in the 
logarithm in the range 
2(R! + R 2 ) to 50AU 
0.4 min (3 min) 
41 min (1400 min) 



Run time (Monte Carlo) 
Run time (NBODY4 with 
GRAPE6Af) 

Notes: where two values are given, the first value refers to runs 
with N a + 2N b = 3750, the second (in brackets) to those with 
N a + 2JV;, = 15000. The timings are on a 3GHz PC (iV-body) and 
AMD Opetron 242 (Monte Carlo). 

energy larger than E tc = —GM/rt are immediately removed 
from the system - M is the total mass and rt is the tidal 
radius. The comparison of the evolution of the mass with 
time (Fig.[T} suggests that a value just about 7 = 0.02 is an 
appropriate choice (especially for an age of order a few Gyr, 
as in M67). 

Apart from mass, the other fundamental measure of a 
cluster is its radius, and the same comparison for the half- 
mass radius (r/J is presented in Fig. [2] It is however, much 
less discriminating of the appropriate value of 7, particu- 
larly for A — 2500. A comparison of the two panels also 
suggests caution in applying the Monte Carlo method to a 
single system with A 10 3 , because of the increasing role 
of statistical fluctuations. 

To properly assess the inferred values of the free param- 
eters of the Monte Carlo code it is important to check the 
intrinsic statistical fluctuation of the code. As can be seen 
from Fig. [3] the spread between models with exactly the 
same parameters, but with different initial random number 
sequence (iseed), is substantial, even for A = 15000. This 
spread is even larger for A = 2500, as can be expected from 
theory. The spread between results with different (3 m in and 
t is well inside the spread connected with different iseed. 
Only the spread between models with different 7 is larger 
that the one connected with different iseed. The best val- 
ues of the free code parameters are: 7 = 0.2, /3 = 0.03 and 
t = 0.01. 



2.2.2 Models with full tidal field 

The process of escape from a cluster in a steady tidal field 
is extremely complicated. Some stars which fulfil the en- 
ergy criterion (binding energy of the star greater than the 
critical energy E tf = -1.5(GM/r t ), see Spitzer (1987)) 
can still be trapped inside the potential well. Those stars 
can be scattered back to lower energy before they escape 
from the system. As was pointed out bv lBaumgardtl (|200ll ) 
these mechanisms cause the cluster lifetime to scale non- 
linear ly with relaxation time, in contrast with what would 
be expected from the standard theory. The efficiency of 
these effects decreases as the number of stars increases. 
To account for this in the Monte Carlo code an additional 
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Figure 3. Comparison of the evolution of the total mass and half- 
mass radius in Af-body and Monte Carlo simulations for models 
with tidal cutoff. The initial conditions are given in Tab.[2land de- 
scribed in the figures. The iV-body model is the heavy continuous 
line, and the others are Monte-Carlo simulations with different 
initial random number sequence iseed = 10, 20, 30, 40, 50. 



free para meter was introduc ed according to the theory pre- 
sented bv lBaumgardtl (|200ll ). The critical energy for escap- 
ing stars was approximated by: Et f = —ot{GM/rt), where 
a = 1.5 - a(ln(-fN)/N) 1/4 . Thus the effective tidal radius 
for Monte Carlo simulations is »"t e ,, = rt/a and it is smaller 
than r t . This leads to the result that for Monte Carlo sim- 
ulations a system is slightly too concentrated compared to 
A-body simulations, but the evolution of the total mass is 
well reproduced, as well as the scaling of the dissolution time 
with A. 

Fig-B]shows the evolution of the total mass and the half- 
mass radius for different a for A = 10000. The other free 
parameters for the case of a full tidal field are the same as 
for the tidal cutoff case: 7 = 0.02, r = 0.01 and fi m in = 0.03. 
As can be seen by comparing Fig[5] (lower two panels) with 
FigC3(top panel), again the spread between models with dif- 
ferent flmin and r is well inside the spread connected with 
different iseed. The statistical spread also does not substan- 
tially interfere with the determination of a and 7 (see FigfS] 
(top panel) for 7). 

As can be seen in Fig. [4] (lower panel) the evolution of 
Th up to time about 0.5 Gyr is slightly too slow in compari- 
son to A-body results. (This effect is even more pronounced 
for smaller N: see Fig[2] top panel.) This behaviour is con- 
nected with the way in which the effect of stellar mass loss 
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(Myr) 



Figure 4. Comparison of the evolution of the total mass and 
the half-mass radius in TV-body and Monte Carlo simulations for 
models with full tidal field. The initial conditions are given in Tab. 
[2] and described in the figures. The TV-body model is the heavy 
continuous line, and the others are Monte-Carlo simulations with 
(from the right) a = 1.0, 1.05, 1.1. 



is fed in to the cluster. In the Monte Carlo model the stellar 
evolution mass loss is postponed until the end of the over- 
all time step, usually several Myr. So, for the most massive 
stars the stellar evolution can be substantially delayed and 
the cluster expands slower. In Fig. [6] the evolution of is 
presented for different models in which the overall time step 
was reduced by factor of two up to a certain time, s. It is 
clear that reduction of the overall time step in the phases 
of cluster evolution in which the most massive stars end 
their evolution helps to bring the Monte Carlo results close 
to the TV-body ones. In later phases of cluster evolution, 
in which the time-scale of stellar evolution becomes larger 
than the half-mass relaxation time, the evolution does not 
depend systematically on the chosen overall time step. (In 
the simulations used in the determination of the free code 
parameters the adopted overall time step was a compromise 
between accuracy and speed.) 

As can be seen in Figs. [7J [8j [9] the Monte Carlo code 
can reproduce TV-body simulations not only in respect of 
the global parameters of the system, but also in respect of 
properties connected with binary activity. Despite the fact 
that the total number of binaries in the system and the 
binary fraction agree quite well with TV-body simulations, 
the total binding energy is substantially too high for the 
Monte Carlo simulations. This is connected with the fact 




t (Myr) 

Figure 5. Comparison of the evolution of the total mass in N- 
body and Monte Carlo simulations for models with full tidal 
field. The initial conditions are given in Tab. [2] and described 
in the figures. The TV-body model is the heavy continuous line, 
and the others are Monte-Carlo simulations with (from the left) 
7 = 0.03,0.02,0.01 (top), (from the ) = 0.06,0.03,0.015 (mid- 
dle), and (from the right) r = 0.02,0.01,0.005 (bottom). 

that the present Monte Carlo simulations cannot follow 3- 
and 4-body interactions directly as the TV-body code does. 
Binaries can only harden or dissolve. Therefore much of the 
complexity of binary dynamical interactions is missing in 
the present Monte Carlo simulations. 

In Fig. [10] is shown the number of collisions as a func- 
tion of time for the Monte Carlo and TV-body models. The 
collisions are mainly connected with binary mergers due to 
stellar evolution or dynamical binary interactions. There are 
only a few direct physical collisions between single stars. Up 
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500 1000 1500 2000 2500 3000 3500 4000 4500 
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Figure 6. Comparison of the evolution of the half-mass radius in 
TV-body and Monte Carlo simulations for models with a full tidal 
field. The initial conditions are given in Tab. [2] and described in 
the figures. The TV-body model is the heavy continuous line, and 
the others are Monte-Carlo simulations with different time, s, up 
to which the overall time step is half of the standard one. The 
model with s = 100 has two different statistical realisations with 
iseed = 20, 40. Note that some curves are overprinted (and hence 
invisible) at early times. 




Figure 7. Comparison of the evolution of the binary binding 
energy in TV-body and Monte Carlo simulations for models with 
full tidal field. The initial conditions are given in Tab. [2] and 
described in the figures. The TV-body model is the heavy con- 
tinuous line, and the others are Monte-Carlo simulations with 
a = 1.1, 1.05, 1.0. 



to time about 1.5 Gyr both models give similar results, and 
then the TV-body model shows a larger number of collisions 
than the Monte Carlo model. Again this can be attributed 
to the fact that in the Monte Carlo simulations the complex 
dynamical binary interactions are not followed. In the Monte 
Carlo code a binary can only coallesce if, after the dynam- 
ical interaction, the periastron distance is smaller than the 
sum of the stellar radii. In the TV-body code binary coalles- 
cence occurs if, during the interaction, the distance between 
two stars is smaller than the sum of their radii. Definitely 
the latter can happen more frequently, in the case of strong 
interactions such as prolonged resonances (temporary cap- 
ture). When most collision events are connected with the 
stellar evolution of nearly isolated binaries, then both mod- 
els agree. 




Figure 8. Comparison of the evolution of the binary fraction in 
TV-body and Monte Carlo simulations for models with full tidal 
field. The initial conditions are given in Tab. [2] and described in 
the figures. The TV-body model is the heavy continuous line, and 
the others are Monte-Carlo simulations with a = 1.1, 1.05, 1.0. 
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Figure 9. Comparison of the evolution of the number of binaries 
in TV-body and Monte Carlo simulations for models with full tidal 
field. The initial conditions are given in Tab. [2] and described in 
the figures. The TV-body model is the heavy continuous line, and 
the others are Monte-Carlo simulations with a = 1.1, 1.05, 1.0. 



2.2.3 Model of M67 

The old open cluster M67 is an ideal testing ground for mod- 
elling the interactions between dynamical and stellar evolu- 
tion. It has a substantial population of blue stragglers, which 
are almost certainly the product of stellar collisions, or merg- 
ers within primordial binaries. Also, it is small enough that 
its enti re life history can be modelled with TV-body tech- 
niques dHurlev et al.|[2005h , though this has become possible 
only within the last few years: though its present mass (of 
order 2OOOM0) makes it seem an easy target for simulation, 
its initial mass is likely to have been much higher (Tab. 
This table also specifies the other initial para meters for our 
model , which closely follow the prescription of Irlurlev et al.l 
(2005). The initial value of the tidal radius is determined by 
scaling the value in their model at 4Gyr to the initial mass 
of our model. 

The da ta from TV-body simulations of M67 

jHurlev et all [2005) were used to calibrate the last re- 
maining parameter of the Monte Carlo code, namely a. 
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Figure 10. Comparison of the evolution of the number of col- 
lisions in TV-body and Monte Carlo simulations for models with 
full tidal field. The initial conditions are given in Tab. [2] and 
described in the figures. The TV-body model is the heavy con- 
tinuous line, and the others are Monte-Carlo simulations with 
a = 1.1, 1.05, 1.0. 
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Figure 11. Mass profiles of the Monte Carlo and TV-body models 
of M67 at 4Gyr, scaled by the total mass involved in constructing 
the profile; from the right all mass, luminous mass (as defined in 
the text) and blue stragglers. The initial conditions are given in 
Tab. [3] and described in the figure. 



Table 3. Initial conditions for M67 



N s + N b 
M(0) 

Initial model 
Initial tidal radius 
Initial mass function 

IMF of binaries 
Binary fraction 
Binary eccentricities 
Binary semi-major axes 

Run time (Monte Carlo) 
Run time (NBODY6) 



24000 

1.904 x 1O 4 M 

Plummer 

32.2pc 

Krouj3a J _Tout_& 



Gilmore (1993 



with ai = 1.3 

Kroupa. Gilmore fc Tout (1991) . eq.(l) 
TV 3 /(TV S + N b ) = 0.5 
/(e) = 2e 

Uniformly distributed in the logarithm 
in the range 2(R 1 + R 2 ) to 50AU 
7 min 
1 month 



The inferred formula is a = 1.5 - 3.0(ln(~/N)/N) 1/4 . 
The comparison of results from TV-body and Monte Carlo 
simulations for M67 confirmed the values of 7, r and f3 m in 
found for smaller TV systems. 

The results of a comparison are summarised in Tab. U 
Taking into account the intrinsic statistical fluctuations of 
both methods the results presented in Tab. [4] show reason- 
ably good agreement. At the time of 4 Gyr, when the com- 
parison was done, both models consist of only a small frac- 
tion of the initial number of stars, making fluctuations even 
more important. The values of the half-mass radius sug- 

Table 4. Monte Carlo and TV-body results for M67 at 4 Gyr 



TV-body (Hurley et al. 20051 This work 



M/Mq 

h 

r t pc 1 
r h pc' 1 
M L /M Q 
M lw /M q 



2037 
0.60 
15.2 
3.80 
1488 
1342 
2.70 



1984 
0.59 
15.1 
3.03 
1219 
1205 
2.67 



L - stars with mass above O.5M0 and burning nuclear fuel 
L10 — the same as L but for stars contained within 10 pc 



gest that the Monte Carlo model is slightly too concentrated 
by comparison with the TV-body model. This, however, can 
be attributed to the treatment of the tide (Sec l2.2.2"j) . which 
leads to a smaller effective tidal radius than the tidal radius 
inferred from TV-body simulations. Additionally, in TV-body 
simulations, stars are considered as escapers on ly if their dis- 
tance from the cluster centre is larger than 2rt l|Hurlev et al.l 
2005). The mass outside rt is about lOOM© (see Fig |ll[) . 
which is small compared to the total cluster mass at any 
time, but nevertheless leads to slightly too large a half-mass 
radius. Compared to TV-body simulations, the values shown 
by the Monte Carlo simulations for the luminous mass Ml 
and the luminous mass inside 10 pc distance from the cen- 
tre Miio are too low. (The luminous mass is the mass of 
all stars with mass above 0.5Mq and burning nuclear fuel 
l|Hurlev et al.ll2005l )). The reasons for this disagreement are 
unclear, but it is partly attributable to the smaller total 
mass of the Monte Carlo model, and balanced by the some- 
what larger total mass in white dwarfs (see below). Note, 
however, that the lower effective tidal radius, r tcff , in the 
Monte Carlo model forces most stars to be confined inside 
10 pc. Therefore Ml and Mlw are practically identical for 
the Monte Carlo simulations. Furthermore, the mismatch 
between the two models is much smaller inside lOpc, which 
suggests that the reason for the disagreement in Ml con- 
cerns mainly large radii. 

The spatial distributions of mass for the Monte Carlo 
and TV-body models are illustrated in Fig. 1111 There are only 
4 blue stragglers in the model, compared with 20 in the TV- 
body model. Agreement is not expected, because our model 
excludes one of the main channels for blue straggler pro- 
duction, i.e. collisions during resonant encounters. Despite 
the large difference in the number of blue stragglers present 
in both models, their mass distributions are very similar. 
Blue stragglers are more centrally concentrated then other 
luminous object in the cluster. The comparison with the TV- 
body model is qualitatively satisfactory, but quantitatively 
reflects the larger half-mass radius and larger effective tidal 
radius in the TV-body model. 

Comparisons of the time-evolution of the number- 
density within the core and half-mass radii for Monte Carlo 
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N-body - Hurley Stars M > 0.8M 0[iot 



3000 
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Figure 12. Number density within the core and half-mass radii 
of the Monte Carlo and TV-body models of M67, as a function of 
time. The initial conditions are given in Tab. |31 and described in 
the figure. 



and TV-body models are presented in Fig. [T^l Agreement 
between the two models is quite satisfactory. The observed 
difference can be attributed to two factors: 

(i) different definitions of the core radius are used in 
the two models. In the Monte Carlo mod el the c o re ra - 
dius is defined according to equation 1-34 in ISpitzer] (Il987l ) 
and in the TV-body mo del from a density-weighted radius 
|Casertano fc Hutlll985l ); 

(ii) the Monte Carlo model is more centrally concentrated 
than the TV-body model because of the smaller effective tidal 
radius r* 



A comparison of the projected number density of lumi- 
nous stars above 0.8Mq for the M onte Carlo and TV-body 
simulations a nd from observation s l|Bica fc Bonattol (120051 ) 
as quoted in Irlurlev et alj (120 05)) is presented in Fig. 1131 
This confirms the conclusions reached so far: the overall den- 
sity of the Monte Carlo model is slightly larger than that of 
the TV-body model and the half-mass radius of the Monte 
Carlo model is too small. Both models exceed the observed 
surface density in the central part of the system but under- 
predict it in the outer halo. These regions require separate 
discussion: 

(i) The higher observational value at large radii sug- 
gests contamination of the observed field by stars which 
ar e not members of the cluster. Indeed, the data quoted 
m iHurlev et all (2005) are not corrected for the background 
density of stars. We also have not corrected the data for the 
background, in order to analyse th e simulation data in the 
same way as m IHurlev et all (|2005l ) (see Fig. 7 there). (We 
consider the effect of the background in Sec[3]|. 

(ii) It is possible that the observational value at small 
radii could be lowered if it were supposed that the obser- 
vations are not fully corrected for the large binary fraction 
in the core. Furthermore, as can be seen from the numer- 
ical simulations, the surface density is lower if only lumi- 
nous stars are taken into account. So, corrections for a low- 
luminosity cutoff and unresolved binaries can help to bring 
observation and simulation closer. However, we should not 
expect a large correction factor for contamination, because 
of the high latitude of M67. 




Figure 13. Profile of the projected number density of all lumi- 
nous stars and binaries and luminous stars above O.8M0, com- 
pared with t h e obs ervations of Bonatto & Bica (extracted from 
IHurlev et alj J2005I . Fig. 7)) and TV-body simulations at 4 Gyr. 
The initial conditions are given in Tab. \3\ and described in the 
figure. 



This discussion of the centre of M67 suggests that probably 
only some changes in the initial model of M67 can bring 
both observations and simulations into agreement, and we 
consider this in Sec[3] (That was not our intention in the 
present section, where our aim is to check that the Monte 
Carlo code can produce results consistent with the TV-body 
model for a realistic cluster model.) 

In Fig. [14] the surface brightness profiles for the Monte 
Carlo and TV-body models are presented. To construct these 
surface brightness profiles all stars and binaries were used. 
The data are very noisy, particularly for the TV-body simula- 
tion. The agreement between the two models is reasonably 
good. Again the conclusion reached before are confirmed. 
The surface brightness in the central parts of the system 
is slightly larger for the Monte Carlo model than that of 
TV-body model and outside in the cluster halo the surface 
brightness is larger for the TV-body model. The latter is again 
connected with the effective tidal radius for the Monte Carlo 
code which is smaller than the nominal tidal radius for the 
two models. Its effects are particularly clear in Fig |13l 

A form of colour-magnitude diagra m is shown in Fig, [T5l 
which can be compared with Fig. 10 of IHurlev et al.l (2005). 
The resemblance is qualitatively satisfactory, except for the 
relative paucity, already referred to, of blue stragglers in the 
Monte Carlo model, and the shortness of the sequence of 
blue stragglers, co mpare d to the TV-body model. 

IHurlev et all (2005) discuss the different exotic popu- 
lations of their model at some length. As already stated in 
connection with the blue straggler population, however, our 
model lacks important processes for the formation of such 
objects. Therefore we confine attention to more normal pop- 
ulations, namely white dwarfs. The mass fraction of white 
dw arfs is 0.18, s l ightly larger than the value of 0.15 given 
by IHurlev et alj l|2005l ). Therefore the total mass in white 
dwarfs is larger in the Monte Carlo model by about 60Mq. 
In particular the white dwarf fraction in the central part of 
the system is larger for the Monte Carlo model than for the 
TV-body model (Fig. I16[l . The spatial distributions of white 
dwarfs are similar in the two models. The maximum lies 
around 6-8 pc, and is more pronounced in the TV-body 
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Figure 14. Surface brightness profile of all luminous stars and 
binaries for the Monte Carlo and TV-body simulations at 4 Gyr. 
The initial conditions are given in Tab. \3\ and described in the 
figure. 



Figure 16. Mass fraction of white dwarfs as a function of radius 
for Monte Carlo and TV-body models at 4 Gyr. The results are 
for all WDs and only single WDs and double WDs. The initial 
conditions are given in Tab. [3] and described in the figure. 



N = 2012 
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B-V 

Figure 15. Colour-magnitude diagram at 4 Gyr for Monte Carlo 
simulation. The initial conditions are given in Tab. [3] and de- 
scribed in the figure. 
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V 

Figure 17. Luminosity function of single main sequence stars for 
Monte Carlo and TV-body models at Gyr and 4 Gyr. The initial 
conditions are given in Tab. [3] and described in the figure. 



model; in the Monte Carlo model the profile could be flat 
(within fluctuations) below this range of radii. The half-mass 
radius of the white dwarfs i s 2.73 pc, much big ger than the 
value of 0.6 pc reported bv lHurlev et al.l l|2005T ); we believe 
their value may be in error. 

For single main sequence stars we present in Fig. I17l the 
luminosity functions for the Monte Carlo and TV-body mod- 
els for times Gyr and 4 Gyr. The luminosity functions for 
time Gyr agree quite well for the two models. Only for 
very bright stars, V < 10 mag, can one observe noticeable 
differences. They can be attributed to statistical fluctuations 
connected with different realisations of the initial model. For 
4 Gyr our Monte Carlo result misses significant numbers of 
stars at the high-luminosity end of the distri bution, which 
are fo und in the TV-body simulation and which iHurlev et"aH 
l|2005l ) attribute to collisions; as already stated in our discus- 
sion of blue stragglers, important channels for the formation 
of such stars are missing at present in our model. For the low 
mass end of the luminosity function the two models agree 
reasonably well. 

Finally, the comparison betw een observations 
^Montgomery. Marschall fc Janes] ll993T ) and the Monte 



Carlo simulations for the luminosity function is presented 
in Fig. [18] The luminosity function was normalised by the 
number of stars with luminosity V < 15.5 mag. In the figure 
are shown the luminosity functions for only main sequence 
stars, only main sequence binaries, and all main sequence 
stars and binaries (for the Monte Carlo model) and for main 
sequence stars and binaries (observations). It is clear that 
the simulations do not agree with the observations. The 
luminosity function for the Monte Carlo model is too low 
at the high-luminosity end, too high at the low-luminosity 
end and too high near the main sequence turn-off (V = 13 
mag). The drop in the luminosity function of the Monte 
Carlo model at the high-luminosity end can be attributed 
to an underproduction of blue stragglers in the Monte 
Carlo model comparable to observations. The excessive 
luminosity function at the low-luminosity end and around 
the main sequence turn off cannot be so easily explained. 
Even though observational issues may be relevant at the 
faint end (Sec(3|, some of these mismatches suggest that 
the initial model of M67 is wrong and some refinement 
is needed. Because the Monte Carlo model agrees so well 
with the TV-body model (Fig |17|) . similar conclusions may 
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Observations Montgomery et al. 1 993 ■ 
MC only MS Slars 

MC only MS Binaries ■ 

MC all MS Stars and Binaries ■ 



M67 

N = 24000 
full tidal field 
f b == 0.5 
m max = 50 M 
|i min a 0.03 
Y= 0.02 

T = 0.01 

«= 1.15 



Table 5. Properties of M67 a 



Figure 18. Comparison of the luminosity function for the Monte 
Carlo model at 4 Gyr and observations (Montgomery at al. 1993), 
for main sequence stars only, binaries only and both main se- 
quence stars and binaries. The initial conditions are given in Tab. 
[3]and described in the figure. Note that all curves are normalised 
separately. 



be reached for a comparison of the TV-body model with 
observed luminosity function, except for the region occupied 
by blue stragglers. 

In this section it has been shown that the Monte Carlo 
code is able to follow the evolution of a realistic star cluster 
model at a similar level of complexity to that of the N- 
body code. The data provided by the code is as detailed 
as for the TV-body code and can be used for comparison 
with many kinds of observational data. The next section 
will be devoted to further Monte Carlo modeling of the old 
open cluster M67 in order to improve the match with its 
observational properties and find possible initial conditions 
for this cluster. 



3 REFINEMENT OF THE M67 MODEL 



As was shown in the previ ous subsection (see Sec. I2.2.3|l the 
model of M67 proposed bv lHurlev etHI ( 2005), whether one 
uses an TV-body code or a Monte Carlo code, shows signifi- 
cant disagreement with such observational properties of the 
cluster as the luminosity function and the surface density 
pr ofile. Other observa tional data about M67 are summarised 
m iHurlev et al] (|2005T ) and listed in Tab. [5] Generally, how- 
ever, the structural cluster parameters are not well known, 
and have been derived from more basic data by some model- 
dependent analysis. Therefore we prefer to compare with the 
observational data as directly as possible. 

The observational data we shall use for comparison 
with the results of Monte Carlo simulations are: (i) the lu - 
minosity function - ([Montgomery. M arschall & Janes 1993), 
and (ii) the surface de nsity profile - (|Bonatto fc BicalboOol ; 
iBica fe Bonattol 120051 ). Unfortunately, the observational 
data for the surface density profile also seem very uncertain, 
even though they were collected by 2MASS (Two Micron All 
Sky Survey); they differ by a factor larger than 1.5. In all fig- 
ures in which the surface density profiles will be presented 
there are two observational curves: (i) - (|Bica fc Bonattol 
l2005h (bottom), corrected for the background density at the 
level of 0.73 stars arcmin -2 which was estimated at the 



Distance from Sun 

Absolute distance modulus b 

Distance from GC 

Orbit eccentricity 

Luminous mass 

Core radius 

Tidal radius 

Half-mass radius c 

Binary fraction 

Z 



870 pc 
9.44 

6.8 - 9.1 kpc 
0.14 

~ 1000M Q 
1.2 pc 
> 11.4 pc 
~ 3.0 pc 
~ 50% 
~ 0.0 
~ 4 Gyr 
0.16 



A v = 3.25E(B - V) 

a References are given in| Hurley et al.l IT2005h 

b (Montgomery, Marschall & Janes 1993); used in the calculation 

of the luminosity function of of the model 

c For main sequence stars with masses 0.87Mq and within 10 
pc 



distance about 25 arcmin, and (ii) - (|Bonatto fc Bicall2003T ) 
(top), corrected for the background density at the level of 4 
stars pc -2 , which was estimated at the limiting radius about 
9 pc. 

The model of M67 presented in the previous section had 
three problems: (i) it was too dense, (ii) it produced too flat 
a luminosity function for dim stars, and (iii) it contained 
too many stars around the main sequence turn off. If we 
take the observations at face value, to bring the model into 
better agreement with observation it has to either lose more 
of its less massive stars, or else contain smaller numbers 
of those stars initially. Additionally, it has to have some 
property such as initial mass segregation or stronger energy 
generation in order to show a smaller concentration at the 
present day. We now describe the parameter space which we 
explored in order to improve the model. 

The free parameters of the initial models are: (i) TV, 
number of objects (stars and binaries), which we explored 
in the range between 22000 and 40000, (ii) fb, binary frac- 
tion, in the range between 0.4 and 0.7, (iii) r±, tidal radius, 
between 30 and 38 pc, (iv) r t /r h , i.e. the ratio between the 
tidal and half-mass radii, between 6 and 12, and (v) cximf, 
the power-law index of the low-mass part of the initial mass 
function (IMF) , between 0.1 and 1.3. (The canonical value 
is olimf = 1.3 (|Kroupall2007l )). Over 135 models of the old 
open cluster M67 were run. We did not carry out this explo- 
ration very systematically; rather on the basis of inspection 
of the results, the parameters of the next set of new models 
were chosen. 

Unfortunately, there is no single model which can re- 
produce the observational properties of the cluster. Instead, 
there are several models which can equally well produce a 
reasonable fit to the observations (Figs. [T^l [20l [2T1 [22l [23|l 
In general from very low values of olimf, i-e. 0.1, up to the 
canonical value, 1.3, the agreement with observations is rea- 
sonably good. For some models the luminosity function is 
modelled better, while for others the surface density profile 
is better. The initial parameters of the best models and the 
parameters at 4 Gyr are summarised in Tab. El The other 
model parameters are close to those chosen bv lHurlev et al.l 
(2005). The free parameters of the Monte Carlo code are 
exactly as determined in the previous Sections. 
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Table 6. Initial parameters, and data at time 4 Gyr, for some models of M67 
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a model with different initial realisation of the sequence of random numbers. All other model parameters are the same as for the model 

above. 
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Figure 19. Comparison of the surface density profile and lumi- 
nosity function for the Monte Carlo model for climf =0.1 and 
observations. The initial model parameters are described in the 
figures. The observational data is described in the text 



Figure 20. Comparison of the surface density profile and lumi- 
nosity function for the Monte Carlo model for ocimf = 0.5 and 
observations. The initial model parameters are described in the 
figures. The observational data is described in the text 



For the purpose of the following comparison with the 
Monte Carlo simulati ons, we shall focus on t he surface den- 
sity profile given in iBica fc Bonatto! (|2005l ). corrected as 
above for the background stars. Clearly, the models show 
reasonably good agreement with the observations (see Figs 
HH]-[22](top panel)). Use of the corrected observational data 
brings the surface density profile for the large radii into much 
better agreement with the simulations (see for comparison 
Fig. II 3p . The special treatment of the tide in the Monte 
Carlo model plays only a minor role: the effective tidal ra- 
dius, rt ef , , is only reduced by about 15% in comparison to 



the true tidal radius, rt. It seems that, despite the large 
latitude of M67, the contamination of the observed surface 
density profile by background stars plays an important role, 
at large radii. 

The comparison between the luminosity functions from 
the Monte Carlo models and observations is given in the 
same figures as for the surface density profiles. The luminos- 
ity functions are in reasonable agreement with observations, 
except that the modelled luminosity function has an excess 
for V > 16 mag, in comparison with observations. Of course 
there are also noticeable differences for the high-luminosity 
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Figure 21. Comparison of the surface density profile and lumi- 
nosity function for the Monte Carlo model for a/MF = 0.7 and 
observations. The initial model parameters are described in the 
figures. The observational data is described in the text 



Figure 22. Comparison of the surface density profile and lumi- 
nosity function for the Monte Carlo model for ciimf — 0.9 and 
observations. The initial model parameters are described in the 
figures. The observational data is described in the text 



end connected with the relative lack of blue stragglers in 
the Monte Carlo models. It could be argued that the ex- 
cess of low-luminosity stars cannot be explained by suppos- 
ing that not all low-mass stars are observed, because the 
observational field of M67 is high-latitude, not heavily con- 
taminated and sparse, and so all stars with V ~ 18 mag 
should be observed. On the other hand the authors them- 
selves (I Montgomery. Marschall fc Janesll 19931 ) declared that 
the determination of the luminosity function "is a difficult 
procedure because of the presence of background stars". 
Their background correction has been applied in the ob- 
servational data shown in these figures, and they did it in 
the following way. The background was estimated by count- 
ing stars in an area of the colour-magnitude diagram just 
blue of the main sequence, and equal in area to the assumed 
boundaries of the main sequence. But photometric errors 
in the colours rise abruptly by V — 15, and the resulting 
spread in the main sequence is very evident below V — 18. 
If the result is that main sequence stars spread into the area 
in which field stars are counted, then the derived (observed) 
luminosity function will be too small. 

Because of the difficulty in quantifying this effect, we 
should consider whether there are dynamical processes, 
omitted from the models, which could also account for the 
mismatch at the faint end of the luminosity function. The 
difference in the luminosity functions for V > 16 mag sug- 
gests some mechanism by which the open cluster M67 very 



efficiently removed a substantial fraction of low mass stars 
(M < 0.7M@) during its evolution. There are two possibili- 
ties: 

(i) - removal of residual gas during the first few million 
years in a cluster with primordial mass segregation. The 
resulting expansion would lead to preferential removal of 
low mass stars. 

(ii) - interactions with the galactic disk and bulge can 
produce tidal shocks which in turn again preferentially re- 
moves low mass stars l|Spitzerlll987M . These effects would, of 
course, alter the entire history of mass-loss in the models, 
and require more massive initial models. 

Finally, we checked the influence of statistical fluctua- 
tions, which are intrinsic to the Monte Carlo method, on 
the determination of the initial cluster parameters which we 
have inferred from the comparison with observations. To as- 
sess the scale of this effect, the same model was repeated 
with different statistical realisations of the initial model, 
(i.e. different initial seeds (iseed) for the sequence of random 
numbers). Typical results are presented in Figs. [2T1 (Model 
III) and EH (Model IV). It is clear that the use of different 
realisations of the initial model has a large impact on the ob- 
servational properties of the cluster at time 4 Gyr. The good 
agreement for the surface density profile within the half- 
mass radius of the cluster is totally destroyed: the new model 
(iseed = 20) has a much higher central surface density, by 
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Figure 23. Comparison of the surface density profile and lumi- 
nosity function for the Monte Carlo model for a/MF = 1-3 and 
observations. The initial model parameters are described in the 
figures. The observational data is described in the text 

factor of about 3. Also the luminosity function for model IV 
shows poorer agreement with the observations than model 
III: there are many more low-mass main sequence stars than 
in model III. As can be seen in the Tab. [6] the model with 
iseed = 20 is less advanced in its dynamical evolution. It has 
larger total mass and half-mass radius, and a smaller binary 
fraction. As can be seen in Fig. 1241 at time 4.5 Gyr the model 
gives similar results to model III. Just because of the differ- 
ent statistical realisation of the initial conditions, model IV 
has smaller average mass than model III, and consequently 
smaller mass loss due to stellar evolution and a correspond- 
ingly smaller cluster expansion. The difference between the 
average masses is only about 4%, but it seems that this is 
enough to change the rate of cluster evolution subst antially. 
This conclusion is supported by Hurley's findings IjHurleyl 
120071 ) that the incidental formation of a massive binary can 
totally change the observational properties of the cluster. 



4 DISCUSSION AND CONCLUSIONS 

In this paper we have presented an advanced Monte Carlo 
code for the evolution of rich star clusters, including most 
aspects of dynamical interactions involving binary and single 
stars, and the internal evolution of single and binary stars. 
It was shown that the free parameters of the Monte Carlo 
code can be successfully calibrated against results of small 




r(po) 



Observations Montgomery et al. 1993 

MC all MS Stars anO Binaries 

t = 4.0 Gyr M c all MS Stars and Binaries at 4.5 Gyr 

M67 




8 10 12 14 16 18 

V 

Figure 24. Comparison of the surface density profile and lumi- 
nosity function for the Monte Carlo model for cijmf — 0.7 and 
observations for different initial realisation of the models, iseed 
= 20. The initial model parameters are described in the figures. 
The observational data is described in the text 



iV-body simulations and simulations of the old open cluster 
M67. 

Sec(3] described our best models of M67. The results 
show that an equally good fit to the observational data can 
be achieved by models which differ substantially in some 
initial parameters, such as the slope of the IMF. By con- 
trast it seems that the other parameters are better con- 
strained, at least within the ensemble of models which we 
studied. Actually, none of the models can successfully fit all 
the observational properties of M67 that we have studied, 
but we have argued that the remaining mismatches can be 
understood in terms of known characteristics of the Monte 
Carlo method, or the observational problem of subtracting 
the background. These difficulties are especially pronounced 
at the bright and faint ends of the luminosity function. The 
most satisfactory models are characterised by the following 
initial parameters: N about 30000, r t about 33 pc, /(, about 
50% and r t /rh about 10. The word "about" is used deliber- 
ately, because of the large effect of statistical fluctuations in 
such small systems. It is worth noting that a satisfactory fit 
can be achieved for a large range of values of the power-law 
index of the IMF for low mass stars, though it seems that 
values of oimf in the range 0.5 - 0.7 give slightly better 
agreement with observations. 

Finally, these Monte Carlo simulations clearly show the 
strong influence of statistical fluctuations on the observa- 
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tional properties of a cluster. In consequence it seems that 
recovery of the initial cluster parameters from a comparison 
between numerical models and observations may be very 
difficult or even impossible, particularly for models with low 
or moderate N. More observational data is needed to con- 
strain the models better, but the data on such properties as 
the number of particular kinds of binaries, pulsars or blue 
stragglers seem to provide only very weak constraints. 

Despite these successes in fitting TV-body models and 
the open cluster M67, the code has some known shortcom- 
ings, which we summarise here. 

(i) Cross sections: it has become apparent 
l|Fregeau fc Rasiol [2007) that the use of cross sections 
can lead to some systematic errors in the evolution of 
core parameters and other quantities. Within the context 
of the Monte Carlo code we are working on, it is known 
how to replace t hese by explicit numerica l calculation of 
the interactions (|Giersz fc Spurzeml [2003), and this will 
be our next improvement. One side effect of the current 
absence of explicit interactions is that we do not yet model 
collisions which occur during them; therefore one channel 
for the formation of blue stragglers is missing from these 
simulations. 

(ii) Higher-order multiples: It is widely argued that pri- 
mordial triples and higher multiples should be incorporated 
into simulations along with primordial binaries. In any case, 
hierarch ical triples form abundantly in binary-binary inter- 
actions (|Mikkolalll984T ). Such higher-order multiples are ig- 
nored in the present Monte-Carlo code, as cross sections for 
interactions with other objects have not yet been devised. 
Hierarchical triples and higher-order multiples can be intro- 
duced as new species when explicit calculation of interac- 
tions has been incorporated. 

(iii) Escape: the Monte Carlo code described here incor- 
porates a tidal cutoff , and a simple mo dification based on 
the theory devised bv lBaumgardtJ (|200ll ). Other treatments 
are possible and worth trying. 

(iv) Rotation: the Monte Carlo code is based on spherical 
symmetry, and would require rather fundamental and very 
difficult reconstruction in order to cope with cluster r ota- 
tion. As was pointed out bv lKim. Lee fc Spurzeml (|2004T ) the 
rotation only somewhat accelerates the rate of core collapse. 

(v) Static tide: t he effect of tidal shocks ha ve been exten- 
sively studied (e.g. iKundic fc Ostriker] l|l995h ) and it would 
be possible to add the effects as another process altering the 
energies and angular momenta of the stars in the simula- 
tions. The addition of tidal shocks will be more important 
when modelling Galactic globular clusters than open clus- 
ters, which usually are confined inside the Galactic disk. 

Despite these limitations, some of which are difficult to 
cure, the Monte Carlo model presented in this paper shows 
its potential power in simulations of star clusters, from open 
clusters to rich globular clusters. Monte Carlo models are 
feasible in a reasonable time (a week or so) for globular 
clusters, which are too large for direct TV-body models, and 
future papers in this series will present results on M4 and 
several other globular clusters. The data provided by Monte 
Carlo simulations are as detailed as those provided by an N- 
body code. No available simulation methods, except Monte 
Carlo and TV-body methods, can really provide that kind of 
comprehensive information. Even when TV-body simulations 



eventually become possible, Monte Carlo models will remain 
as a quicker way of exploring the parameter space for the 
large scale TV-body simulations. 
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